The ARMA model is a stationary model, which means we will need to assume (and verify) that our data are a realization of a stationary process. Because regression can be used to break down a nonstationary model to a trend, seasonal components, and a residual series, it’s often reasonable to apply stationary models to the residuals from a time series regression.
A time series is called a mixed autoregressive moving average model of order \((p,q)\) if it is stationary and takes the form:
\[ x_t = \phi_1 x_{t-1} + ... + \phi_p x_{t-p} + w_t + \theta_1 w_{t-1} + ... + \theta_q w_{t-q}\] We can observe that the model includes a white noise component, an AR component, and a moving average component. This equation assumes the series has already been demeaned.
We can rewrite this equation in terms of backshift operators and simplify in order to get a very concise version of the equation.
\[ x_t - \phi_1 x_{t-1} - ... - \phi_p x_{t-p} = w_t + \theta_1 w_{t-1} + ... + \theta_q w_{t-q}\] \[ x_t (1-\phi_1 B - ... - \phi_p B^p) = w_t (1 + \theta_1 B + ... + \theta_q B^q)\] \[ x_t (1-\phi_1 B - ... - \phi_p B^p) = w_t (1 + \theta_1 B + ... + \theta_q B^q)\]
\[ \phi_p(B) x_t = \theta_q(B) w_t\]
Important points about an \(ARMA(p,q)\) model:
If we wanted to accommodate a non-zero mean \(\mu\), then we would add a constant term on the right-hand side: \(\alpha = \mu (1-\phi_1 - ...- \phi_p)\):
\[ x_t = \alpha + w_t +[ \phi_1 x_{t-1} + ... + \phi_p x_{t-p} ] + [\theta_1 w_{t-1} + ... + \theta_q w_{t-q}]\]
The equation for \(\alpha\) can also be rearranged to show that the mean \(\mu\) is stationary:
\[ \mu = \frac{\alpha}{1-\phi_1 - ...- \phi_p}\]
If we go back to the ARMA equation in terms of the backshift operator:
\[ \phi_p(B) x_t = \theta_q(B) w_t\]
we see that we can rearrange it to get an MA process (\(x_t\) as a function of \(w_t\)):
\[ x_t = \frac {\theta_q(B)}{ \phi_p(B)} w_t \equiv \psi(B) w_t \] Through expansion of \(\phi_p(B)\) in the denominator as a geometric series, this becomes an \(MA(\infty)\) model.
Likewise, we can rearrange to achieve an \(AR(\infty)\) model (\(w_t\) as a function of \(x_t\)):
\[ w_t = \frac{ \phi_p(B)} {\theta_q(B)} x_t \equiv \pi(B) x_t \]
Skipped. Terminology in async is completely inconsistent from slide to slide. It seems like it’s mainly used to derive an equation corresponding to the form of the ACF plot.
Reminder about the general properties of the ACF and PACF plots:
| Plot | AR(p) Model | MA(q) Model | ARMA (p,q) Model |
|---|---|---|---|
| ACF | Tails off | Abrupt cutoff after lag q | Tails off |
| PACF | Abrup.t cutoff after lag p | Tails off | Tails off |
Because the ACF for \(ARMA(1,1)\) and \(AR(1)\) only differ by a constant, it’s hard to use ACF plots to distinguish between the models. In an \(AR(p)\) model, as \(p\) approaches 1, the series gets more persistent. However, even with a smaller \(p\), the addition of an \(MA(q)\) component will also contribute to persistence.
We may be able to more easily distinguish (in some cases) from a PACF plot. Recall that the PACF plot for a \(AR(p)\) model cuts off abruptly after \(p\) lags; a PACF plot for an \(ARMA(p,q)\) plot generally will take longer to tail off.
In practice, we will estimate several models and use AIC, BIC, or forecasting a test set to decide between them.
df <- read.table('pounds_nz.dat', header=TRUE)
bpnz <- ts(df$xrate)
str(bpnz)
## Time-Series [1:39] from 1 to 39: 2.92 2.94 3.17 3.25 3.35 ...
head(bpnz)
## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## [1] 2.924 2.942 3.172 3.254 3.348 3.507
tail(bpnz)
## Time Series:
## Start = 34
## End = 39
## Frequency = 1
## [1] 3.050 3.187 3.229 3.192 3.352 3.531
par(mfrow=c(2,2))
plot(bpnz)
hist(bpnz)
acf(bpnz)
pacf(bpnz)
Observations:
We’ll try an MA model anyways as a demo:
(ma5 <- arima(bpnz, order=c(0,0,5))) # MA(5) model
##
## Call:
## arima(x = bpnz, order = c(0, 0, 5))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 intercept
## 1.539 1.334 1.169 0.645 0.27 2.873
## s.e. 0.196 0.296 0.277 0.265 0.17 0.112
##
## sigma^2 estimated as 0.015: log likelihood = 24.17, aic = -34.34
AIC(ma5)
## [1] -34.34
BIC(ma5)
## [1] -22.7
The first 4 coefficients are statistically significant, but we still need to do diagnostics:
ma5r <- resid(ma5)
par(mfrow=c(2,2))
plot(ma5r)
hist(ma5r)
# qqnorm(ma5r)
acf(ma5r)
pacf(ma5r)
Observations:
We also do the Ljung-Box test for autocorrelation of the residual series. The null hypothesis is of independence (no correlation):
Box.test(ma5r, type="Ljung-Box") # using default of 1 lag
##
## Box-Ljung test
##
## data: ma5r
## X-squared = 0.5, df = 1, p-value = 0.5
The high p-value says we cannot reject the null hypothesis of independence.
Now let’s look at in-sample vs. out-of-sample fit:
par(mfrow=c(1,1))
ts.plot(bpnz, fitted(ma5), resid(ma5), lwd=c(1,2,1),
lty=c('solid','dashed','dashed'), col=c('black','blue','black'),
main='MA(5) model')
The in-sample fit looks reasonable. Here’s what an \(MA(5)\) model would forecast:
# do a forecast of 6 steps
(fcast <- forecast(ma5, h=6))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 40 3.507 3.348 3.666 3.263 3.750
## 41 3.396 3.107 3.684 2.954 3.837
## 42 3.239 2.882 3.596 2.693 3.785
## 43 3.051 2.650 3.453 2.438 3.665
## 44 2.921 2.507 3.335 2.288 3.554
## 45 2.873 2.457 3.289 2.237 3.509
plot(fcast)
lines(fitted(ma5), lty='dashed', col='blue')
Let’s try back-testing for out-of-sample performance:
bpnz_bt <- window(bpnz, end=33) # hold out last 5 observations for testing
ma5_bt <- arima(bpnz_bt, order=c(0,0,5))
fcast_bt <- forecast(ma5_bt, h=12)
par(mfrow=c(1,1))
plot(fcast_bt)
lines(window(bpnz, start=34), lty='solid', col='black')
lines(fitted(ma5_bt), lty='dashed', col='blue')
Observations:
Next we’ll try an ARMA model for the sake of demo. (The requirement of stationarity is not actually met for this time series.)
(arma11 <- arima(bpnz, order=c(1,0,1))) # ARMA(1,1) model
##
## Call:
## arima(x = bpnz, order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.892 0.532 2.960
## s.e. 0.076 0.202 0.244
##
## sigma^2 estimated as 0.0151: log likelihood = 25.14, aic = -42.27
AIC(arma11)
## [1] -42.27
BIC(arma11)
## [1] -35.62
arma11r <- resid(arma11)
par(mfrow=c(2,2))
plot(arma11r)
hist(arma11r)
acf(arma11r)
pacf(arma11r)
Box.test(arma11r, type="Ljung-Box") # using default of 1 lag
##
## Box-Ljung test
##
## data: arma11r
## X-squared = 0.014, df = 1, p-value = 0.9
Observations:
This looks like a better model for this data.
Now let’s look at in-sample fit and forecast.
par(mfrow=c(1,1))
ts.plot(bpnz, fitted(arma11), resid(arma11), lwd=c(1,2,1),
lty=c('solid','dashed','dashed'), col=c('black','blue','black'),
main='ARMA(1,1) model')
In-sample fit looks reasonable, although generally shifted slightly to the right of actual.
(fcast <- forecast(arma11, h=6))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 40 3.532 3.375 3.690 3.292 3.773
## 41 3.471 3.197 3.744 3.052 3.889
## 42 3.416 3.077 3.755 2.898 3.934
## 43 3.367 2.984 3.750 2.781 3.952
## 44 3.323 2.908 3.738 2.689 3.957
## 45 3.284 2.846 3.722 2.614 3.954
plot(fcast)
lines(fitted(arma11), lty='dashed', col='blue')
Observations:
Now, back-test, holding back 6 observations:
arma11_bt <- arima(bpnz_bt, order=c(1,0,1))
fcast_bt <- forecast(arma11_bt, h=12)
par(mfrow=c(1,1))
plot(fcast_bt)
lines(window(bpnz, start=34), lty='solid', col='black')
lines(fitted(arma11_bt), lty='dashed', col='blue')
Observations:
ARIMA models are one way of dealing with nonstationary time series, which may have trends or seasonal effects. Simple “differencing” can off convert a nonstationary time series to a stationary one. Often a first-order differencing is sufficient (but generally we don’t want to go higher than 2nd order).
Note: not all nonstationary time series can be deal with by differencing. Notably, volatility clustering (conditional heteroskedasticity) that’s common in financial times series requires a different kind of model. Those are commonly model with Autoregressive Conditional Heteroskedastic (ARCH) models.
The term “integrated” arises from the fact that a differenced series needs to be aggregated in order to recover the original. The simplest \(I(0)\) process is white noise; the simplest \(I(1)\) process is the random walk (because after first differencing, we have a white noise model).
For example, a random walk has the following form:
\[ x_t = x_{t-1} + w_t\]
Rearranged, we get a stationary white noise series \(w_t \sim N(0,\sigma_w^2)\):
\[ x_t - x_{t-1} \equiv \nabla x_t = w_t \]
Just a quick refresher as to random walk time series look like. Note that the drift term plays the same role as the slope in deterministic linear trend models.
par(mfrow=c(2,1))
# without drift
x <- w <- rnorm(100)
for (t in 2:length(x)) {
x[t] <- x[t-1] + w[t]
}
plot(ts(x), main='Random Walk Simulation')
# with drift
x <- w <- rnorm(100)
del <- 0.5
for (t in 2:length(x)) {
x[t] <- del + x[t-1] + w[t]
}
plot(ts(x), main='Random Walk Simulation w/ Drift (delta=0.5)')
As a reminder, the expectation value grows over time due to the drift:
\[E(x_t)=x_0 + t \delta\]
and the variance grows without bounds:
\[ Var(x_t) = t \sigma^2\]
The autocovariance simplifies to:
\[\gamma_k = t \sigma^2\]
Because autocovariance is a function of time, this model is obviously nonstationary.
First differencing can also remove deterministic trends. If we have linear trend of the following form, we can consider either first differencing (which results in an \(MA(1)\) progress) or simply subtracting the trend (and analyzing residuals):
\[ x_t = a + bt + w_t\] \[ \nabla x_t = x_t - x_{t-1} = (a - a) + (bt - b(t-1)) + (w_t - w_{t-1}) = b + w_t - w_{t-1}\] Reducing and rewriting in terms of the lag operator, we have a stationary \(MA(1)\) process:
\[ \nabla x_t = b + \theta_q(B) w_t\] Alternately, we could have subtracted the trend (instead of differencing) to achieve a white noise process:
\[ x_t - (a + bt) = a + bt + w_t - (a + bt) = w_t\] If our original time series showed increasing variance over time, we should try a log transformation (and then differencing, if there’s also a trend we want to try to remove).
An \(ARIMA(p,d,q)\) model is an \(ARMA(p,q)\) model that’s applied after taking the \(d^{th}\) difference of the original time series \(x_t\). Using the lag operator, this can be expressed as:
\[ \phi_p(B)(1-B)^d x_t = \theta_q(B) w_t\]
Let’s simulate this model:
\[ x_t = 0.5 x_{t-1} + x_{t-1} - 0.5 x_{t-2} + w_t +0.3 w_{t-1}\] Rearranging: \[ x_t - x_{t-1} = 0.5(x_{t-1} - x_{t-2}) + w_t + 0.3 w_{t-1}\] \[ x_t - x_{t-1} - 0.5(x_{t-1} - x_{t-2}) = w_t + 0.3 w_{t-1}\] \[ \nabla x_t - 0.5 \nabla x_{t-1} = w_t - w_{t-1}\] Since we have \(x_t\) and \(x_{t-1}\) terms, we can rewrite in terms of the lag operator. This has the form of an \(ARIMA(1,1,1)\) model:
\[ \nabla x_t (1 - 0.5B) = (1 + 0.3B)w_t\] \[ \nabla x_t = 0.5B \nabla x_t + (1+0.3B)w_t\] \[ \nabla x_t - 0.5 \nabla x_{t-1} = w_t + 0.3 w_{t-1} \] We see that after applying the first differencing operator \(\nabla\), the \(ARIMA(1,1,1)\) model is transformed to an \(ARMA(1,1)\) stationary model with an AR parameter \(\phi_1 = 0.5\) and an MA parameter of \(\theta_q=0.3\). (Be careful of signs!)
Let’s simulate this equation, and then compare plots between the original and the differenced equation:
x <- w <- rnorm(100)
for (t in 3:100) {
x[t] = 0.5 * x[t-1] + x[t-1] - 0.5 * x[t-2] + w[t] + 0.3 * w[t-1]
}
x_ts <- ts(x)
# alternately, we could have asked arima.sim to run the simulation
# x_ts <- arima.sim(model=list(order=c(1,1,1), ar=0.5, ma=0.3), n=100)
par(mfrow=c(1,2))
plot(x_ts, main='Original Simulated Time Series')
plot(diff(x_ts), main='Differenced Simulated Time Series')
acf(x_ts, main='Original Simulated Time Series')
acf(diff(x_ts), main='Differenced Simulated Time Series')
pacf(x_ts, main='Original Simulated Time Series')
pacf(diff(x_ts), main='Differenced Simulated Time Series')
Observations:
Now let’s estimate the model and plot the original and fitted time series.
(arima111 <- arima(x_ts, order=c(1,1,1)))
##
## Call:
## arima(x = x_ts, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## 0.476 0.248
## s.e. 0.125 0.130
##
## sigma^2 estimated as 0.912: log likelihood = -136.2, aic = 278.4
(arma11_diffed <- arima(diff(x_ts), order=c(1,0,1)))
##
## Call:
## arima(x = diff(x_ts), order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.475 0.248 0.044
## s.e. 0.126 0.130 0.226
##
## sigma^2 estimated as 0.912: log likelihood = -136.2, aic = 280.4
par(mfrow=c(1,1))
ts.plot(x_ts, fitted(arima111), lwd=c(1,2),
lty=c('solid','dashed'), col=c('black','blue'),
main='ARIMA(1,1,1) model')
We see essentially the same results when we manually difference before fitting an ARMA model. Our estimated coefficients are quite close to the true (supposedly unknown) process that we simulated. Our fitted values from the model are quite close to the originals.
Now let’s do our standard diagnostics with the model’s residuals:
arima111r <- resid(arima111)
par(mfrow=c(2,2))
plot(arima111r, main='Residual TS')
hist(arima111r)
acf(arima111r)
pacf(arima111r)
Box.test(arima111r, type='Ljung-Box')
##
## Box-Ljung test
##
## data: arima111r
## X-squared = 0.014, df = 1, p-value = 0.9
The plots and Ljung-Box test both fail to reject independence (ie. supporting the residual TS as a realization of a white noise process), indicating that we should be able to forecast with this model.
Let’s attempt forecasting:
(fcast <- forecast(arima111, h=12))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 101 6.316 5.0922 7.540 4.4442 8.188
## 102 6.051 3.6118 8.491 2.3204 9.782
## 103 5.925 2.4087 9.441 0.5473 11.303
## 104 5.865 1.4110 10.319 -0.9467 12.677
## 105 5.836 0.5596 11.113 -2.2337 13.906
## 106 5.823 -0.1865 11.832 -3.3675 15.013
## 107 5.816 -0.8548 12.487 -4.3862 16.019
## 108 5.813 -1.4639 13.090 -5.3161 16.942
## 109 5.812 -2.0265 13.650 -6.1757 17.799
## 110 5.811 -2.5515 14.173 -6.9783 18.600
## 111 5.811 -3.0456 14.667 -7.7337 19.355
## 112 5.810 -3.5135 15.134 -8.4493 20.070
par(mfrow=c(1,1))
plot(fcast)
lines(fitted(arima111), lty='dashed', col='blue')
We observe that while the model fit is good, the forecast is essentially a flat line.
These are basically an extension of ARIMA models that add a lag equal to a number of seasons in order to remove seasonal effects. The form of the model is \(ARIMA(p,d,q)(P,D,Q)_m\) where \(p,d,q\) are the nonseasonal lag terms, and \(P,D,Q\) are the seasonal lag terms, where \(m\) is the number of periods per year (e.g. \(m=12\) for monthly seasonality, or \(m=4\) for quarterly seasonality). The \(P,D,Q\) terms basically just represent backshifts of a seasonal period. We use lowercase \(\phi_p,\theta_q\) to represent nonseasonal components, and \(\Phi_P, \Theta_Q\) to represent the seasonal terms:
\[ \Phi_P(B^m) (1-B^m)^D \phi_p(B) (1-B)^d x_t = \Theta_Q (B^m) \theta_q(B) w_t\] These terms correspond to:
\[ [\text{Seasonal AR(P)}][\text{Seasonal Diff}][\text{Non-seasonal AR(p)}][\text{Non-seasonal Diff}] = [\text{Seasonal MA(Q)}][\text{Non-seasonal MA(q)}]\]
For, example a monthly \(ARIMA(1,1,1)(1,1,1)_4\) model looks like:
\[ (1-\Phi_1 B^4) (1-B^4)^1 (1-\phi_1 B) (1-B)x_t = (1+ \Theta_1B^4)(1+\theta_1B)w_t\] ### SARIMA Examples
This model indicates that a monthly value is affected by the same month’s value in the previous year, is:
\[x_t = \alpha x_{t-12} + w_t\]
The characteristic equation for that formula is:
\[ (1 - \alpha B^{12}) = 0\] Rearrange, we get:
\[ B = \frac{1}{\alpha}^{1/12} = \alpha^{-1/12}\] The model is stationary when \(|\alpha^{-1/12}|>1\).
\[ x_t = x_{t-1} + \alpha x_{t-12} - \alpha x_{t-13} + w_t\] We could also write as: \[ \nabla x_t = \alpha \nabla x_{t-12} + w_t\] which helps us understand the intuition that the change at time $ t $ depends on the change at the same time of the previous year. Rearranging and factorizing gives:
\[ \Theta_1(B^{12})(1-B)x_t = (1-\alpha B^{12})(1-B)x_t = w_t\] We see the seasonal AR(1) term, with \(m=12\), and the nonseasonal difference term. This model is nonstationary, since we see that there is a root of \(B=1\).
A quarterly, seasonal moving average model \(ARIMA(0,0,0)(0,0,1)_4\) (stationary, without a trend) is:
\[ x_t = (1-\beta B^4)w_t = w_t - \beta w_{t-4}\] If a nonseasonal stochastic trend was also present, a \(ARIMA(0,1,0)(0,0,1)_4\) model (using 1st order difference to remove the trend) could look like:
\[ x_t = x_{t-1} + w_t - \beta w_{t-4}\] If, instead, a seasonal stochastic trend was present, the \(ARIMA(0,0,0)(0,1,1)_4\) model could be used to remove the seasonal stochastic trend:
\[ x_t = x_{t-4} + w_t - \beta w_{t-4}\] ### Worked example: EU Retail Trade Index (Quarterly)
library(fpp) # samples for Rob Hyndman's book
## Loading required package: fma
## Loading required package: expsmooth
## Loading required package: lmtest
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: tseries
## Warning: package 'tseries' was built under R version 3.4.4
data(euretail)
str(euretail)
## Time-Series [1:64] from 1996 to 2012: 89.1 89.5 89.9 90.1 89.2 ...
Plot the time series:
par(mfrow=c(3,1))
# original time series
plot(euretail)
# take seasonal difference (still not stationary)
plot(diff(euretail, lag=4))
# take nonseasonal and seasonal differences (better)
plot(diff(diff(euretail, lag=4)))
Now let’s try modeling it. Note that if we use sarima(...), then we need a slightly different form to retrieve the residuals. The function comes with a nice, default display.
library(astsa)
##
## Attaching package: 'astsa'
## The following object is masked from 'package:fpp':
##
## oil
## The following objects are masked from 'package:fma':
##
## chicken, sales
## The following object is masked from 'package:forecast':
##
## gas
eu_arima <- sarima(euretail, 0,1,1,0,1,1,4)
## initial value -0.629077
## iter 2 value -0.801933
## iter 3 value -0.814397
## iter 4 value -0.827539
## iter 5 value -0.829110
## iter 6 value -0.829179
## iter 7 value -0.829179
## iter 8 value -0.829179
## iter 8 value -0.829179
## iter 8 value -0.829179
## final value -0.829179
## converged
## initial value -0.830649
## iter 2 value -0.831162
## iter 3 value -0.831170
## iter 4 value -0.831171
## iter 4 value -0.831171
## iter 4 value -0.831171
## final value -0.831171
## converged
eu_arima$ttable # display estimated parameters
## Estimate SE t.value p.value
## ma1 0.2901 0.1118 2.595 0.012
## sma1 -0.6909 0.1197 -5.771 0.000
eu_arima$AICc # AIC, BIC are also available
## [1] -0.608
par(mfrow=c(2,2))
eu_arimar <- resid(eu_arima$fit)
plot(eu_arimar)
hist(eu_arimar)
acf(eu_arimar)
pacf(eu_arimar)
Since both the ACF and PACF show significant spikes at lag \(k=2\), this suggests we need to include additional nonseasonal terms. (We don’t see signifiant seasonal lags.) Also observe that the Ljung-Box plot (nice!) rejects the null hypthesis of independence of the residuals. So, we need to do more work on a better model. Next, I’ll go straight to the best model that the book found:
#eu_arima <- sarima(euretail, 0,1,3,0,1,1,4)
# eu_arima$ttable # display estimated parameters
# eu_arima$AICc # AIC, BIC are also available
# eu_arimar <- resid(eu_arima$fit)
(eu_arima <- arima(euretail, order=c(0,1,3), seasonal=list(order=c(0,1,1), period=4)))
##
## Call:
## arima(x = euretail, order = c(0, 1, 3), seasonal = list(order = c(0, 1, 1),
## period = 4))
##
## Coefficients:
## ma1 ma2 ma3 sma1
## 0.263 0.370 0.419 -0.661
## s.e. 0.124 0.126 0.130 0.156
##
## sigma^2 estimated as 0.145: log likelihood = -28.7, aic = 67.4
par(mfrow=c(2,2))
eu_arimar <- resid(eu_arima)
plot(eu_arimar)
hist(eu_arimar)
acf(eu_arimar)
pacf(eu_arimar)
Box.test(eu_arimar, type='Ljung-Box')
##
## Box-Ljung test
##
## data: eu_arimar
## X-squared = 0.003, df = 1, p-value = 1
NOTE: I used arima instead of sarima because the latter produced an object that I couldn’t pass to forecast.
Observations:
Now, let’s try a forecast:
(fcast <- forecast(eu_arima, h=12))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2012 Q1 95.21 94.72 95.70 94.47 95.96
## 2012 Q2 95.29 94.50 96.07 94.08 96.49
## 2012 Q3 95.38 94.26 96.50 93.67 97.09
## 2012 Q4 95.40 93.90 96.91 93.11 97.70
## 2013 Q1 94.63 92.73 96.53 91.72 97.54
## 2013 Q2 94.64 92.39 96.90 91.20 98.09
## 2013 Q3 94.64 92.06 97.23 90.69 98.60
## 2013 Q4 94.67 91.75 97.58 90.21 99.12
## 2014 Q1 93.89 90.61 97.17 88.87 98.91
## 2014 Q2 93.91 90.28 97.53 88.36 99.46
## 2014 Q3 93.91 89.94 97.88 87.84 99.98
## 2014 Q4 93.93 89.62 98.24 87.34 100.52
par(mfrow=c(1,1))
plot(fcast)
lines(fitted(eu_arima), lty='dashed', col='blue')
Observations:
First load and take a quick look at the dataset.
cbe <- read.table('cbe.dat', header=TRUE)
head(cbe)
elec <- ts(cbe$elec, start=1958, freq=12)
str(elec)
## Time-Series [1:396] from 1958 to 1991: 1497 1463 1648 1595 1777 1824 1994 1835 1787 1699 ...
head(elec)
## Jan Feb Mar Apr May Jun
## 1958 1497 1463 1648 1595 1777 1824
tail(elec)
## Jul Aug Sep Oct Nov Dec
## 1990 14002 14338 12867 12761 12449 12658
summary(elec)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1463 3239 5891 6312 8820 14338
quantile(elec, c(0.01,0.05,0.1,0.25,0.5,0.75,0.9,0.95,0.99))
## 1% 5% 10% 25% 50% 75% 90% 95% 99%
## 1597 1848 2108 3239 5891 8820 11332 12192 13844
Now look at the plots.
par(mfrow=c(2,2))
plot(elec)
hist(elec)
acf(elec)
pacf(elec)
The first difference of removes the trend, but not the variance. First difference of log transformed time series removes both:
par(mfrow=c(3,1))
plot(elec)
plot(diff(elec), main='First Difference of TS')
plot(diff(log(elec)), main='First Difference of Log(TS)')
We can see the improvement quantitatively:
summary(elec)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1463 3239 5891 6312 8820 14338
summary(diff(elec))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1471.0 -197.0 -30.0 28.3 258.5 1383.0
summary(diff(log(elec)))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.14533 -0.03832 -0.00675 0.00540 0.05105 0.18836
Now look at the ACF and PACF of the differenced and log-transformed dataset:
elec_trans <- diff(log(elec))
par(mfrow=c(1,2))
acf(elec_trans, lag.max=48)
pacf(elec_trans, lag.max=48)
Both still show seasonal effects, made even more clear by extending the lag axis. I can’t duplicate the results in the async slides, so I’ll instead replicate the best model found in the book. (See Cowper page 144 for an example of a function to determine the best parameters.)
(elec_arima <- arima(log(elec), order=c(0,1,1),
seasonal=list(order=c(2,0,2), period=frequency(log(elec)))))
##
## Call:
## arima(x = log(elec), order = c(0, 1, 1), seasonal = list(order = c(2, 0, 2),
## period = frequency(log(elec))))
##
## Coefficients:
## ma1 sar1 sar2 sma1 sma2
## -0.640 0.481 0.514 -0.079 -0.471
## s.e. 0.046 0.238 0.236 0.223 0.142
##
## sigma^2 estimated as 0.000417: log likelihood = 956.8, aic = -1902
elec_arimar <- resid(elec_arima)
par(mfrow=c(2,2))
plot(elec_arimar)
hist(elec_arimar)
acf(elec_arimar)
pacf(elec_arimar)
Box.test(elec_arimar, type='Ljung-Box')
##
## Box-Ljung test
##
## data: elec_arimar
## X-squared = 0.017, df = 1, p-value = 0.9
Observations:
Now, try forecasting:
(fcast <- forecast(elec_arima, h=12))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1991 9.442 9.416 9.468 9.402 9.482
## Feb 1991 9.412 9.384 9.440 9.369 9.454
## Mar 1991 9.492 9.462 9.521 9.447 9.537
## Apr 1991 9.433 9.402 9.464 9.386 9.480
## May 1991 9.532 9.499 9.564 9.482 9.581
## Jun 1991 9.574 9.540 9.608 9.523 9.625
## Jul 1991 9.610 9.575 9.645 9.557 9.664
## Aug 1991 9.606 9.570 9.642 9.550 9.661
## Sep 1991 9.511 9.474 9.548 9.454 9.568
## Oct 1991 9.511 9.472 9.549 9.452 9.570
## Nov 1991 9.478 9.439 9.518 9.418 9.539
## Dec 1991 9.488 9.448 9.529 9.426 9.551
par(mfrow=c(1,1))
plot(fcast)
lines(fitted(elec_arima), lty='dashed', col='blue')
Observations:
We need to manually exponentiate fitted and forecasted values if we them on original scale:
# exponentiate the forecasted values
exp_fitted <- exp(fcast$fitted)
exp_fcast <- exp(fcast$mean)
exp_upper95 <- exp(fcast$upper[,'95%'])
exp_lower95 <- exp(fcast$lower[,'95%'])
exp_upper80 <- exp(fcast$upper[,'80%'])
exp_lower80 <- exp(fcast$lower[,'80%'])
# ts.plot(elec, exp_fcast, exp_fitted, exp_lower95, exp_upper95, lwd=c(1,2,1,1,1),
# lty=c('solid','dashed', 'dashed','solid','solid'), col=c('black','blue','blue','gray','gray'),
# main='Exponentiated Forecast (Original Units)', xlim=c(1980,1992))
# manually recreate the forecast plot (with CI) using my exponentiated series
plot(elec, xlim=c(1980,1992), ylim=c(6000,16000), main='Exponentiated Forecast (Original Units)')
tmp <- seq(from = 1991, by = 1/12, length = 12) # need x axis for drawing polygon
polygon(x= c(tmp,rev(tmp)),y= c(exp_upper95,rev(exp_lower95)), col="lightgray", border=NA)
polygon(x= c(tmp,rev(tmp)),y= c(exp_upper80,rev(exp_lower80)), col="lightsteelblue3", border=NA)
lines(exp_fcast, lwd=2, lty='solid', col='blue')
# lines(exp_upper95, lwd=1, lty='dotted', col='blue')
# lines(exp_lower95, lwd=1, lty='dotted', col='blue')